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Abstract. Cardiac diseases are the principal cause of human morbidity 
and mortality in the western world. The electric potential of the heart is a 
highly complex signal emerging as a result of nontrivial flow conduction, 
hierarchical structuring and multiple regulation mechanisms. Its proper 
accurate analysis becomes of crucial importance in order to detect and 
treat arrhythmias or other abnormal dynamics that could lead to life- 
threatening conditions. To achieve this, advanced nonlinear processing 
methods are needed: one example here is the case of recent advances in 
the Microcanonical Multiscale Formalism. The aim of the present paper 
is to recapitulate those advances and extend the analyses performed, spe- 
cially looking at the case of atrial fibrillation. We show that both ECG 
and intracardial potential signals can be described in a model-free way as 
a fast dynamics combined with a slow dynamics. Sharp differences in the 
key parameters of the fast dynamics appear in different regimes of tran- 
sition between atrial fibrillation and healthy cases. Therefore, this type 
of analysis could be used for automated early warning, also in the treat- 
ment of atrial fibrillation particularly to guide radiofrequency ablation 
procedures. 
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1 Introduction 

Heartbeat signals emerge as the result of a nontrivial and highly complex syn- 
chronization process between the cardiac pacemaker cells which are hierarchi- 
cally connected. This structure ensures the robustness of heartbeats and under 
normal conditions it exhibits a regular main rate perturbed by small chaotic 
fluctuations. Typically, the amplitude of such fluctuations is much smaller than 
the average interbeat interval, something that makes the healthy (sinus rhythm) 
heartbeat appear as mainly periodic. Nevertheless, the fluctuations around this 
main period are not an unstructured random noise but follow a complex dynam- 
ics. Moreover, extensive research findings in the last decade have shown that, 
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the characterization of these fluctuations can be crucial for determining whether 
the heart is healthy or it is indicating signs of a transition to an arrhythmia, 
despite still expressing normally |16|17|24|28I38) . 

The cardiac action potential is leaded by polarization of pacemaker cells. 
These cells are not homogeneously distributed but mainly concentrate in nodes 
(sinoatrial and atrioventricular nodes), bundles and the Purkinje fibers which 
innervate all the ventricular myocardium. The action of pacemaker cells controls 
the cardiac contractions (atrial and ventricular systole) and relaxation (dias- 
tole) in an organized way to ensure the optimal pumping 22 . Polarization and 
depolarization of the membrane potential is a collective effect and so they are 
affected by the conduction dynamics of the electric flow (2D]- I n this sense, the 
orientation of cardiomyocites |41) plays an important role and other structural 
particularities such as that of epicardial fat |4I5) can have an important effect in 
particular cases. The complex mechanical structure of the heart together with 
the conductivity inhomogeneities and the hierarchical activation are themselves 
governed through nontrivial regulation mechanisms, which include the sympa- 
thetic and parasympathetic divisions of the autonomic nervous system and the 
endocrine system. 

All these elements and, most importantly, their mutual interrelations make 
the cardiac dynamics highly complex and as a consequence its description re- 
mains mainly phenomenological and the microscopic state of the components is 
impossible to determine noninvasively in most cases. That is why effective de- 
scriptions provide an alternative for the automated analysis of heartbeat signals. 
In that context, first studies of interbeat fluctuations found that they have a mul- 
tiscale structure |21l23j and so fractal models were proposed for it [31] . Later on, 
the development of more advanced analysis techniques based on multiresolution 
analysis and characterization of singularities, large deviations and predictability 
permitted a more extensive study [16], which shows that a healthy heartbeat 
has a multifractal structure, while a heart under congestive heart failure devi- 
ates from multifractality. In that study, the multifractal analysis performed is 
based solely on the Legendre spectrum. As a consequence, from a thermodynamic 
point of view, it corresponds to a canonical formalism, not microcanonical, and 
so local singularity exponents for each point cannot be directly obtained within 
the formalism. For that reason, recent methods developed under the framework 
of the Microcanonical Multiscale Formalism (MMF) 35 32 52] make them an 
especially appropriate alternative for heartbeat analysis; e.g. in [38137) . 

The process resulting in the multiscale structure observed in heartbeat con- 
sists of the synchronization process in a hierarchic complex network [40]; the 
connectivity network of cardiac pacemaker cells in our case. As a consequence, 
an analysis based on the singularity exponents [26|25j and the optimal wavelet 
[36 34] applied to heartbeat time series allows directly accessing the geometric 
features that characterize their multiscale behavior. The results obtained in [To] 
are based on a canonical analysis, meaning that the behavior of statistical av- 
erages is used to indirectly retrieve the geometric features: scaling exponents of 
partition functions estimate a curve that can be used to obtain the so-called 
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singularity spectrum by means of a numerically-estimated Legendre transform) . 
This methodology is known to give less accurate estimation on the tails of the 
singularity spectrum for which a microcanonical analysis has been found to be 
much more robust and accurate |50] . Having such estimation has a capital im- 
portance for anticipating as much as possible when heartbeat dynamics starts 
drifting from the healthy behavior. Given the quickness with which heart failure 
can be fatal or leave irreversible sequelae, the precise estimation provided by the 
MMF has a strong potential in helping to save lives and improve the health of 
people with cardiac diseases. Multifractal models originated from the study of 
turbulent flows. While blood can be in turbulent regime -and notably inside the 
heart- it is unlikely that this turbulence is reflected to the cardiac electric flow. 
This electric flow, and in consequence the measured heartbeat signal, are the 
result of a complex synchronization process between pacemakers. As explained 
in [3J, traffic in a complex network under certain regimes becomes multifrac- 
tal and so the analysis of the pacemaker network could explain the observed 
mult if r act ality. 

This paper is presented as an extended version of the preliminary results by 
the same signing authors that were presented in the 6th International Conference 
MDA 2011 in New York 38 . Here we summarize our previous findings, we 
expand the technical description of the methods, we include the presented results 
and we extend them with a deeper analysis, a more exhaustive data processing 
and statistical testing. 

The paper is structured as follows: in the next Section we present the methods 
used in our analysis. We introduce the basics of the MMF and the algorithms 
to accurately retrieve the empirical singularity exponents from a signal. We 
also present a reconstruction method and how it can be used to separate the 
fast dynamics implied by the singularity exponents from a slow dynamics that 
indicates changes in regime. The fast dynamics is a simple orientation transition 
without memory. The slow dynamics modulates the fast one and it can be easily 
sifted from it. In Section [3J we introduce the empirical data to be analyzed as 
well as some basics of the atrial fibrillation condition that are relevant for its 
signal processing. In Section [4] we apply that analysis to the heartbeat data and 
discuss how it can be used to identify dynamical changes, specially for the case 
of atrial fibrillation. Finally, in Section [5] we draw the conclusions of our work. 

2 Methods: The Microcanonical Multiscale Formalism 

The MMF is a theoretical and methodological framework for the analysis of 
multiscale signals. Its basic element of description is by means of the singularity 
exponents of a signal under analysis, which are the exponents describing the 
local regular/singular behavior of the signal around each point [35132152] . 

2.1 Singularity exponents 

Singularity exponents have different mathematical definitions depending on the 
context they are used. The usual notion in complex-signal analysis is related 
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to the Holder or Hurst exponents, including their respective generalizations. 
Although different definitions are possible, the conceptual goal is always the 
same: to describe how the function evolves around a given point by converging 
to a value (regular) or diverging (singular). 

In the most general case, given a signal s that is defined on R d domain and 
images to a R m space, the Holder exponent h(x) of point x is the exponent 
satisfying the following limit, when it exists [18] : 

|| s(x + r) - s(x) || = a(x)r ft W + ( r M*)) ( r ) (1) 

where r = ||r||. This means that in the proximity of x the signal follows a power 
law of exponent ft.(x). An alternative definition that analytically is slightly more 
restrictive is usually called the Hurst exponent |44ll9j and defined as s(x + r) — 
s(x) = (a(x)|r) r-^O)- 1 -). ( r ' l ( x )) where a(x) is a continuous (1, 1) tensor. For 
the purpose of this article, analysis of ID signals of 1 component the definitions 
actually coincide. 

The concept of singularity exponent can be interpreted also in terms of dif- 
ferentiability. A function that is strictly n-derivable at point x has a singularity 
exponent /i(x) = n. So that in this sense the singularity exponent can be related 
to non-integer differentiability. In a similar way, as we will see below, it is also 
related to the content of information. 

Nevertheless, Holder or Hurst exponents defined this way have very specific 
applicability (e.g., in the case of multiaffinc functions) and cannot be directly 
found in real-world signals. The main reason is that the basic power-law be- 
havior is masked by the presence of long-range correlations, noisy fluctuations, 
discretization and finite-size effects. All these make that the analytical limit 
described is not practically attainable |47|52j . and a generalized definition of 
singularity exponent is needed. To achieve this, the objective is to find a certain 
measure /i for which we could take a similar limit: 

A*(S r (x)) = a(x) r d+h{x) + o (V +,l(x) ) (r^O) (2) 

where d is the dimension of the domain, i.e., d = 1 in the ID case, and B r {~x) is 
a ball centered around x having a radius r for a certain norm (choice to be done 
for multi-dimensional cases; they all coincide in ID). 

The actual definition of singularity exponent that we will be using in the 
present article works well in practice and is little affected by the artifacts men- 
tioned above. For it, we will work on the gradient-modulus measure of the signal 
[4"T] . This measure is defined from its density: 

d/i(x) = ||V S ||(x)dx (3) 

a definition that is absolutely continuous with respect to the Lebesgue measure. 
Hence, the measure of any Borelian A is given by: 



p{A) = f dx||V s ||(x) 

J A 



(4) 
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The gradient-modulus measure characterizes the local singularity of any 
point. A signal that has a Holder exponent h(x) + 1 according to eq. (JT]) will 
fulfill also eq. ([2]), with this +1 shift. 

Practical calculations of eq. ^ can benefit from using wavelet-projected 
interpolations, this way effectively avoiding some of the discretization effects 
[6 27; . The wavelet projection of the measure at point x and scale r is expressed as 
7s?/z(x, r) = L d d/x(x') r~ d & ((x — x') /r) with & being a predetermined function 
known as the mother wavelet. As we can see, the operator Iq, is a map from the 
set M. of cx-finite measures on R d to the set of functions M d x IR + — > M. That is 
why a signal that has a singularity exponent at the point x according to eq. {2} 
exhibits this same exponent when wavelet-projected [6 47 , i.e, 

7^//(x,r) = Q*(x) r h(x > + o(> (x) ) (r -> 0) (5) 

It is worth mentioning that wavelet projections expressed in this way treat 
the wavelet function as a kernel for the measure and no additional restriction 
is imposed. This way, we are not limited to use only admissible wavelets (i.e., 
wavelets that form a basis of a function space). In particular, we can use always- 
positive kernels that do not have zero-crossings. High-order wavelets that exhibit 
several zero-crossings have a significant loss in spatial resolution |46|50j . but 
positive kernels minimize spatial spread and can normally reach the original 
resolution, that is, one pixel in the original signal. 

2.2 Reconstruction formula 

The existence of a multiscale hierarchy details dynamical redundancies in the 
signals that can be exploited to reconstruct them from partial information |51j . 
In the following, we describe a reconstruction algorithm for multiscale signals 
that is model-agnostic and only assumes empirically verifiable hypotheses on 
them. 

The starting point is equation ([5]) : the signals under analysis exhibit a singu- 
larity exponent h(x) at each point [3019] . In addition, the values of the singularity 
exponents are organized so that they define a hierarchy of multiscale geometri- 
cal structures |S]- In the literature, these two hypotheses are commonly verified 
statistically at a global level [8143] but they can be expressed also at each point 
if using a specially appropriate multiscale functional [51147152] . The reconstruc- 
tion formula aims at reconstructing the complex signal from the vertex of the 
singularity exponent hierarchy. 

That reconstruction is proved from the theoretical point of view in the case 
of multifractal functions when the vertex of the hierarchy is associated to the 
so-called Most Singular Component (MSC) [43], which is the set comprising 
the points with most singular (i.e., most negative) values of h(x). Therefore it 
is reasonable to study its reconstruction capabilities on empirical signals, es- 
pecially when that component constitutes a meaningful attractor of the signal 
dynamics. This can be achieved in practice by analyzing the information content 
in the signal to consider the most informative set of points as its MSC. Provided 
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that the reconstruction works, the most informative set actually constitutes an 
Unpredictable Point Manifold [39]. Signal predictability is a core concept in the 
processing of complex signals |2|lj , where measures such as Lyapunov exponents 
and Kolmogorov-Sinai entropy are commonly used to identify information con- 
tent. 

The reconstruction algorithm is designed to act on the gradient of the signal 
(or the derivative in the case of ID functions we present here) because the 
multiscale functional used, eq. ([3]), is gradient-based and so this point simplifies 
the development. In order to keep the algorithm as general as possible, we look for 
a universal operator. This operator reconstructs a signal from its Unpredictable 
Point Manifold - ideally, from its MSC [51] - and it is consistent with the known 
statistical invariances of multiscale signals [9] . 

The algorithm is model-agnostic in the sense that it does not assume any 
particular distribution or correlation structure for the MSC set of points or 
the associated singularity exponents. The only assumptions on the projecting 
operator is that it is required to be deterministic, linear, translational invariant, 
isotropic and leading to the known power-spectrum shape. For the case of time 
series, there is an alternative description that substitutes the isotropy condition 
by the causality condition. These requirements uniquely define the operator |51j . 

For a given multiscale signal s let us denote by J 7 ^ its MSC, i.e., the set of 
points x such that /i(x) €=]ft.oo — A, + A[. is the minimum value of h(x) 
over the signal and A is a small, infinitesimal threshold (ideally). The essential 
gradient of the signal is defined as V s(x) = Vs(x) 5^ (x). In that context, 5^ 

means a delta distribution over the J-^ (it uniformly weighs MSC points and 
vanishes outside of the MSC). 

The reconstruction formula introduced in [51j can be expressed as: 

S (x) = (ff-V a)(x) (6) 

where ( • ) is the convolution dot-product and the vector field g is the recon- 
struction kernel. That kernel can be expressed in Fourier space in a compact 
way: 

= ipja (7) 

with i = y— 1 as the imaginary unit, k as the frequency vector and the hat * 
indicating Fourier transform. 

Furthermore, the reconstruction formula determines how the MSC is related 
to the most informative set of points in the signal. The derivation in |51j implies 
that any set T reconstructing the signal must verify: 

div (V s) = (8) 

where T c is the complementary set oiJ-. That applies in particular to the MSC, 
but since the divergence operator is local and eq. (|8]) is linear, the reconstructibil- 
ity of a point can be decided based on its neighborhood only. Therefore, the MSC 
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points are those that must be included to J- to have a full reconstruction of the 
signal, because their values cannot be predicted just knowing the values in their 
surroundings. In other words, the MSC points are the unpredictable points of 
the signal [52] . 

3 Atrial fibrillation 

Atrial fibrillation ( AF) , the most common form of cardiac arrhythmia, is respon- 
sible for significant morbidity each year in all parts of the world. It results from 
the chaotic operation of the top of the heart (atria). Although a priori it is not 
a severe condition by itself, it causes a high mortality rate by its most severe 
complications, either from heart failure or by stroke-related embolism. The main 
cause of AF is related to a change in the electrical conduction properties in some 
type of cardiac tissue: some areas of this tissue depolarize spontaneously or slow 
the spread of the pulse. Today the primary treatment for AF remains medica- 
tion, but ineffective drugs, intolerance to them or their side effects has led to the 
development of new forms of treatment, mainly radiofrequency ablation. This 
ablation technique requires the introduction of a catheter within the heart in 
order to burn the areas that have become electrically deficient. A low- voltage 
radiofrequency current is applied to the pathogenic areas heating them to induce 
the necrosis of the tissue at tiny width and depth (in the order of few mm). In 
cases of paroxysmal AF, Hai'ssaguerre et al. have shown [12] that for 80% of 
patients, the pathogenic tissue is located in one of the 4 pulmonary veins, and 
electrical insulation, obtained by surgical means, allows the patient to regain a 
normal heart rhythm |7|42|45|14|15|13|TT] . But in more severe cases, i.e., per- 
sistent or permanent AF, locating the pathogenic areas remains difficult and is 
still an unsolved problem. During the ablation procedure, catheters inserted in 
the heart can analyze finely the electrical activity of the atria. The morphology 
of the signals obtained and their temporal evolution must guide the surgeon to 
the location of the sources to ablate. But the complexity of the acquired signal 
makes analysis very difficult and there is no clear identification of important 
information. It is in this context where nonlinear analysis techniques like the 
MMF can be applied to identify the changes in cardiac regime that lead to the 
recovery of the normal sinus rhythm. 

3.1 The data 

The processed data consist of multilead recordings of the electric potential mea- 
sured on six human patients affected by different types of atrial fibrillation who 
have undergone radiofrequency ablation in their atrial endocardia. Measures 
have been taken before, during and after RF ablation at the Haut-Leveque hos- 
pital in Pessac. There are both unipolar and bipolar leads, and there are ECG 
measurements on the skin concurrently with electrode catheters inside the heart. 
The standard ECG leads recorded are I, II and III (bipolar) and VI (unipo- 
lar); these constitute the first 4 data channels in our datafiles. Three electrode 
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catheters are introduced through the main veins of the patient. One of them is 
the radiofrequency catheter, which is used for the ablation but is able to measure 
the potential as well; it consists of four electrodes and so provides potential mea- 
sures at two positions, namely the distal one and the proximal one. The second 
one is called pentaray and it is a catheter that can be opened at its end dividing 
in five branches that spread equiangular ly; each branch has four electrodes mea- 
suring at a distal and a proximal positions. These two catheters are introduced 
inside the heart. A third catheter containing ten electrodes is introduced to the 
coronary sinus and provides measures at five positions. In total we have 4 ECG 
potential measures on the skin and 17 intracardial measures for each time in- 
stant. Signals are sampled at a rate of 1 kHz for a total of 1 183 232 data points. 
Figure Q] shows radiographics illustrating the electrode placements. 




Fig. 1. Radiographic images taken during the radiofrequency ablation procedure and 
showing the three catheters and the electrode placement. In both cases, we can observe 
an electrode catheter measuring along inside the coronary sinus, the ablation catheter 
and the pentaray catheter taking measures. The ablation catheter is placed curved and 
was not operating at the moment of the radiographies (it is nevertheless always taking 
measures also). The pentaray catheter is shown taking measures in the left inferior 
pulmonary vein area (left image) and in that of the right superior pulmonary vein 
(right image). 

For illustration purposes, we show in Figure [2] how the fibrillation (desyn- 
chronized beat) is seen in the right atrium and how AF is still slightly manifested 
in the measure from the electrode on the skin VI. 

As a complement, we have also processed data from MIT-BIH Arrhythmia 
Database [29110) . These consist of long and extensive ECG measurements under 
different arrhythmic regimes, in many cases containing also the transition period 
between regimes. The most remarkable particularity is probably the fact that 



9 



0.0 



-0.5 



-1.0 





Fig. 2. Illustration of normal and atrial fibrillation heartbeats. At the top, we show 
external electrocardiogram measures, namely from VI electrode, for the case of sinus 
rhythm (left) and atrial fibrillation (right). At the bottom, we show the measurements 
of the first pair of electrodes of the catheter on the endocardial area around the right 
superior pulmonary vein. Fibrillation is clearly seen inside the heart, and can be still 
noticed outside. 
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they are annotated at the beat level by multiple experts independently (and 
these annotations have been audited and verified by the full community of the 
database users since 1980). They consist of two channels: a bipolar lead (II) and 
a unipolar V* lead, sampled at 360 Hz for 30 min. for a total of 47 patients 
of Boston's Beth Israel Hospital between 1975 and 1979. Electric signals were 
originally recorded analogically into magnetic tapes, but are today distributed 
as digital files. We have processed files 200 to 205 as those appeared to be the 
ones with most AF episodes. 

The significant differences between one dataset and the other made us expect 
important differences in their analyses, particularly in case spurious effects had 
appeared. On the contrary, we have observed consistent results, which indicates 
the robustness and meaningfulness of the analyzing methodology performed. 

4 Results: Heartbeat analysis 

The first step to perform on the empirical signals consists on the validation of 
the mult ifr act ality hypothesis under the MMF. Previous studies have validated 
multifractal dynamics under a canonical formalism, which is based on statistical 
averages (the scaling of structure functions or partition functions). However, the 
MMF is directly geometrical instead of statistical, which means that we charac- 
terize the multiscale character at each point of the signal. In Figure [3] we show 
the resulting singularity spectrum estimated for the ensemble of signals |33l38j . 
Convergence to the presented curve has been observed for all the individual 
signals between the empirical error bars. 

Once the heartbeat data has been validated to be microcanonical, i.e., it 
fits the singularity-exponent relation seen in eq. ([5]), the next step consists in 
analyzing the dynamical properties and characterizing transition points. A first 
observation is that the measures inside the heart clearly have more abundance of 
most singular points. This is something that we expected is because the difference 
observed - more singular points - means that the internal signals are more 
informative and require more points to be reconstructed. Multifractal signals can 
be reconstructed from their component of most singular points (which therefore 
concentrates all the information) , if such component exists |43j . In our case, the 
most singular component (i.e., the most informative) is also larger and more 
dispersed in the internal than in the external measures. We show below that 
the studied signals are effectively reconstructible though nevertheless, the key 
dynamical parameters detected in internal measures are still noticeable in the 
external measures. 

A general theory and methods about reconstructibility of multifractal sig- 
nals from the Most Singular Component (MSC) can be found in [52] ■ With the 
MMF we have access to the singularity value of each point, which is the basis 
of the multiscale coordination of the components in the signal. Additionally, the 
methodological tools provided by this formalism are effectively adapted to work 
with real signals (with discretization, noise, artifacts, aliasing, correlations), and 
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Fig. 3. Singularity spectrum showing the Hausdorff dimension of exponent level sets in 
front of singularity exponent values for the processed heartbeat data. As expressed in 
Section [2 the computed singularity exponents refer to the gradient-modulus measure 
and so singular (negative) values do not imply discontinuities in the signal but in its 
gradient, as there are no exponents below -1. In addition, the full-support component, 
i.e. that of dimension 1, corresponds to a positive singularity value, thus ensuring the 
integrability of the trajectory. 
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we can use this information to look at the dynamical properties of the heartbeat 
series. 

4.1 A dynamical model for analysis 

Given the points we have, we propose a dynamical model for the analysis of 
heartbeat signals. The first key element comes from the calculation of the singu- 
larity exponents and, in particular, their associated Most Singular Component 
(MSC). An interesting aspect of the MSC is that we do not need the values of the 
signal over its points nor the actual values of singularity exponent. In fact, most 
of the dynamical information of the signal is contained only in the orientation 
over the MSC [48] and we can use an analysis of this orientation as the starting 
point. 

First of all, we calculate the singularity exponents of the signals to determine 
where their MSC is located. To achieve this, we locate where are the smallest 
(most singular) exponents in the series, giving a small tolerance level to account 
the numerical fluctuations. Then we define the oriented MSC as a function that 
is zero everywhere except for the most singular points, where it takes the value 
of the sign of the gradient. 

When we apply the reconstruction formula defined in |51j to the oriented 
MSC, the result is a reduced signal that coincides with the original signal at 
short scales. At long scales, there is a slow divergence between them. This allows 
to model the dynamics as a combination of a fast dynamics driven by the MSC 
orientation and a slow- varying field that acts as a factor on it [49j . 

More concretely, we define the oriented MSC as 5oo(*) taking +1 on MSC 
points of positive derivative, -1 on MSC points of negative derivative and 
on non-MSC points. This way, transitions from one point to the other can be 
described as a Markov chain. This way, we call a the Markov states. The two- 
point joint probability is noted as: 

F(a ,a T ) = (P(S oo (t) = a ,S oo (t + T) = a T )) (9) 
and as a consequence the marginal probabilities are: 

PM = (P(M*) = °*))t = (p(^°t))^ (io) 

Under the hypothesis of distributional stationarity, we can expand the process 
as transitions between all the states [48] . The transition at two steps is expressed 
as: 

P(a 2 \a ) = J2 p (^Wi)P(a 1 \a ) (11) 

Cl 

i.e., twice the one-step transition. This means that we can represent all the 
Markovian dynamics through P(a\\ao). The matricial expression of the process 
gives the so-called transition matrix. In our case, with three possible states: +1, 
-1,0. 
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00 0+ 

+0 ++ + 
-0 -+ - 



) 



(12) 



In a Markov process, a state is called recurrent when the probability of start- 
ing at that state and returning to it after a finite number of steps is one. A 
stationary distribution for the process exists if and only if all the states are re- 
current and all the expected times of first return are finite. This is always the case 
in our process, because all transitions are possible and none is absorbing. That 
stationary marginal distribution is a state of convergence of the process. Apply- 
ing T to it results the same, so this is the first eigenvector of the process with 
an eigenvalue of 1. The other two secondary eigenvalues give a characterization 
of the dynamics. 

Hands on the data: The first calculation to be done is the validation of eq. 
(TTTj) : for all the signals the differences between the empirical two-steps P(o2\vo) 
couple values and those derived from combination of two P(<7i|<7o) steps have 
a median relative discrepancy of 1%. These differences are compatible with the 
empirical estimation sampling errorbar and mean that in all cases the memory of 
the MSC orientation decays very fast so that the Markovian behavior hypothesis 
for it could be reasonably assumed. 

A preliminary look at the oriented MSC distributions does not show sharp 
differences between a signal and another. The estimation of the eigenvalues of 
T matrices are primarily limited by their least probable elements. Sampling 
error propagation is around 7% for each lead and heartbeat regime in the Haut- 
Leveque database and 3% for each lead and heartbeat regime in the MIT-BIH 
arrhythmia database. 

To this effect, we have analyzed the Markov processes for all the data, classi- 
fied in four categories: internal channels under Atrial Fibrillation (AF) internal 
channels under (normal) sinus rhythm, ECG channels under AF and ECG chan- 
nels under sinus rhythm. We did not observe significant differences from one 
patient to the other or from one channel to the other inside the category, so 
we have grouped them to enhance presentation of the results and maximize the 
precision. That said, when processing one single channel the results are already 
stable, so we conclude that the method is robust and little data demanding. 

As we can see in Figure 0] the particular signature of AF is conserved when 
the signal inside the heart is propagated to the skin. We notice that despite the 
signals clearly differ in appearance, the dynamical parameters in terms of MSC 
orientation dynamics for the internal measures are still externally observed with 
minimal differences, while at the same time the AF condition significantly per- 
turbs these values both internally and externally. The same behavior is confirmed 
at processing the MIT-BIH arrhythmia database signals. This suggests that it 
should be possible to finely monitor the AF evolution and severity from external 
ECG measures by means of advanced statistical measures that are stable and 
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robust. Even more, transitions to and from fibrillation could be immediately 
detected. 



Signal group 


Second eigenvalue 


Second- eigenvalue error 


Third eigenvalue 


Third-eigenvalue error 


AF skin 


0.28 


0.04 


0.25 


0.04 


AF intracardial 


0.28 


0.02 


0.18 


0.02 


AF MIT-BIH 


0.29 


0.03 


0.19 


0.03 


N skin 


0.19 


0.04 


0.16 


0.04 


N intracardial 


0.19 


0.02 


0.14 


0.02 


N MIT-BIH 


0.18 


0.03 


0.14 


0.03 



Table 1. Values for 2nd and 3rd eigenvalues and their respective errors as plotted in 



Figure [3] 



Source field: As we have mentioned above, the oriented MSC drives the fast 
dynamics, meaning that the sign clustering determines a big part of the dynam- 
ical structure. Nevertheless, on a longer scale the MSC alone does not precisely 
reconstruct the signal because the constant factor is not really constant but 
slowly evolves. This leads to the definition of the source field. Given a signal 
s and the reduced r constructed from its oriented MSC, the source field p is 
defined such that: 

Vs(.i) = p(x) Vr(x) (13) 

but nevertheless this definition is not usable in practice in points outside the 
MSC, because the reduced signal is not well defined in them. To solve this, we 
make use of a more rigorous definition of source field in terms of measures [49] : 

fi s (A) = / dp r {x)p{x) (14) 

J A 

which implies that the source field p(x) is the Radon-Nikodym derivative of the 
two measures: p{x) = dp s /dp r . 

There exist several different strategies to numerically estimate a Radon- 
Nikodym derivative. Since in our case we are most interested in the detection of 
slow dynamical transitions, we have used an iterative algorithm that fits eq. (|14[) 
in a piecewise constant fashion. This way, we concentrate on the determination 
of the dynamical borders and let the MSC lead all the fast Markovianly-stablc 
evolution. The results, shown in Figure [5] show that the source field varies infre- 
quently though it exhibits quite sharp transitions. In consonance with the results 
for the MSC orientation in the subsection above, the dynamical character in the 
case of AF is significantly different than under sinus rhythm. 

We observe a correspondence of the transition points with the points in which 
the deviation between the original and the reconstructed series is more impor- 
tant. So these transition points correspond to transitions in the reconstructibility 
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Fig. 4. Transition matrices T of the oriented MSC have been calculated for each signal 
in each regime; we show their respective second and third eigenvalues which are the 
parameters leading that orientation dynamics for the atrial fibrillation cases (top panel) 
and the normal sinus rhythm (middle panel) . Correspondences between channels show 
that we can group all the internal and all the external leads for each regime thus further 
reducing the sampling error; that is shown in the bottom panel with 6 signal groups: 
Haut-Leveque internal (intracardial) measures, Haut-Leveque ECG measures and MIT 
signals, each for healthy and AF cases. Numerical values are given in Table [T] 
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Fig. 5. Time evolution of the source fields and the reconstructions of an illustration 
signal corresponding to measures from VI electrode under atrial fibrillation. At the 
left, source field displayed (solid) over the original signal (dashed). At the right, signal 
reconstruction (dashed) based on the source field and the Markov-chain modeled MSC. 
Signs of Atrial Fibrillation are noticed in the dynamical parameters. Reconstruction is 
of pretty good quality, especially for the peaks. 



and in the content of information, which means that the detected transitions cor- 
respond to actual changes in the dynamical properties of the signal. The concrete 
mechanism that establishes the link of correspondence with some electrophysio- 
logical transition is nontrivial and complex. 



5 Conclusions 



We have performed a detailed study of cardiac fluctuation signals under the 
approach of the Microcanonical Multiscale Formalism (MMF). The structural 
complexity inherently present in those signals challenges their analysis. The 
emergence of a complex dynamic behavior hinders most attempts to carefully 
estimate microscopic dynamical parameters in the system, especially unless car- 
diac models are assumed with heavy calibration or data assimilation. Therefore, 
an analysis describing in terms of effective dynamics is especially appropriate 
for heartbeat. These signals have been known to be of multiscale character since 
several decades ago [21I23I31] but surprisingly enough this aspect is presented 
most of the times in the literature as coincidences between separate scale levels, 
without exploiting the capabilities of multiscale structure and mutual relation- 
ships - of course with notable exceptions such as in [16]. The description in 
terms of multiscale effective dynamics is essentially model-free, and it assumes 
hypotheses on the data that are easily verifiable. In that sense, this approach is 
close to that of Machine Learning methodologies, but has the advantage that it 
tracks closely the key exponents of the emergence of complexity in the system 
and so the parameters with a physical meaning, be it in terms of information 
content, information transfer, dynamical attractors or critical sets |33j . 
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The multiscale structure in heartbeat is reflected as a definite geometrical 
distribution around manifolds of singularity. This way, the signal under analy- 
sis is decomposed into different components depending on their characteristic 
singularity exponent. This fact decomposes the signal into separate regimes ac- 
cording to their characteristic dynamics. The value of the singularity exponent 
characterizes the power-law behavior under scale changes and directly indicates 
the information content of the component. Consequently, our analysis provides 
direct access to the dynamical structure at each point of the signal and it unveils 
the geometry of singularity components and characterizes the degree of infor- 
mation contained in them. When further exploited, this analysis shows that the 
most singular component (MSC) contains the information of the entire signal 
and can restore it. In other words, this component drives the dynamics of the 
signal. 

A first observation is that we reproduce under the MMF the same type of 
multifractality characterizations for heartbeat series that have been reported 
under a canonical framework [16] . Moreover, the expected observation that data 
from catheters inside the heart have a richer, more informative, more singular 
multiscale signature than those taken on the skin, does not seem to affect the 
MSC parameters of the signal or its reconstructibility, which means that the key 
dynamical features can still be detected on the skin. 

In addition, we have shown that the characteristic dynamical parameters 
retrieved from MMF analysis -namely the orientation of the MSC- can be dy- 
namically described as a Markov chain. Furthermore, these dynamic parameters 
(in particular, the transition eigenvalues) under atrial fibrillation (AF) are sig- 
nificantly different from the sinus-rhythm case. Additionally, they can be equally 
detected from intracardial electrodes or from standard electrocardiogram mea- 
sures on the skin. Therefore, a possible application would be early detection of 
transitions to or from the AF. The main database processed consists of potential 
measures taken during AF ablation procedures at Haut-Leveque Hospital, and 
the robustness of the approach is confirmed by similar results obtained for the 
MIT-BIH Arrhythmia Database which contains significantly different signals. 

Finally, we have noticed that the signal reconstructibility implies that it 
can be separated into the fast MSC-orientation process and a complex slowly- 
varying source field that modulates it. These source fields accurately describe 
the multifractal dynamic changes, what would suggests a possible relationship 
with transitions in electrophysiological processes, such as the evolution of the 
cardiac regulatory mechanism and changes in conductivity of the tissue, as well 
as other structural changes and drifts. 
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